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Introduction 


Over  the  past  few  years,  the  authors  have  been  engaged  in  a 
coordinated  research  program  on  inverse  problems.  By  this  time,  the 
results  have  been  spread  over  a number  of  different  papers.  Thus,  it  is 
the  considered  opinion  of  the  authors  that  a review  of  the  present  state 
of  this  research  be  written.  This  review  also  provides  an  opportunity 
to  filter  some  of  the  work  that  has  proven  to  be  of  less  interest  and  to 
incorporate  into  the  presentation  new  insights  that  have  been  developed 
during  this  research  project. 

This  report  describes  the  three  main  segments  in  this  research 

program: 

(i)  physical  optics  far  field  inverse  scattering  (denoted  by  the 
acronym,  POFFIS); 

(ii)  seismic  or  subsurface  profiling  in  media  with  small  variations 
in  propagation  speed; 

(iii)  analysis  of  the  inverse  source  prob''  m. 
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1.  The  POFFIS  Identity 

The  POFFIS  identity  was  originally  derived  by  BOJARSKI  [1]. 
Analysis  of  the  basic  identity  and  its  extensions  can  be  found  in  LEWIS 
[2],  BOJARSKI  [3],  PERRY  [4],  TABBARA  [5,6],  ROSENBAUM-RAZ  [7],  BLEISTEIN 
[8,9]  and  MAGER  and  BLEISTEIN  [10].  The  problem  of  interest  here  is 
reconstruction  of  the  image  of  a convex  scatterer  from  observations  of 
the  high  frequency  far  field  scattering  from  the  object  in  response  to  a 
known  probing  signal. 

The  POFFIS  identity  relates  the  scattered  field  to  the  Fourier 
transform  of  the  characteristic  function  of  the  scatterer.  The  character- 
istic function  is  equal  to  unity  in  the  region  occupied  by  the  scatterer 
and  zero  elsewhere.  Thus,  knowledge  of  this  function  makes  possible  the 
reconstruction  of  an  image  of  the  scatterer.  The  Fourier  transform  of 
the  characteristic  function  is  itself  a function  of  a vector  transform 
variable,  J<.  The  magnitude  of  the  vector  is  proportional  to  frequency; 

I the  direction  of  the  vector  is  determined  by  the  source-receiver  directions. 

I Consequently,  in  any  practical  situation,  with  limited  bandwidth  and 

directions  of  observation,  the  Fourier  transform  is  known  only  in  some 
aperture  limited--band  limited  and  aspect  angle  1 imi ted— domain  in  J< 
space.  The  problem  of  extracting  information  about  the  characteristic 
I function  from  a high  frequency  aperture  limited  Fourier  transform  is 

. solely  a question  in  the  Fourier  analysis  of  piecewise  constant  (in  this 

i case  zero-one)  functions.  The  problem  is  discussed  from  this  point  of 

I 

|j  j|  view  in  Section  1.3  and  will  prove  useful  in  the  discussions  in  Sections 


2,1— 2.3  below. 
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1.1  Derivation  of  the  POFFIS  Identity 

A point  source  located  at  x in  Figure  1,  gives  rise  to  a 

• 0 

signal,  Ut(x,  x , oj),  which  is  a solution  of  the  reduced  wave  equation, 

i “ “0 

2 2 2 

(1.1)  7 Ur  + Ui  uJc  = - 6(x  - X );  x = (x  , x , x ). 

^ ' I 1 “ -0  " 1 2 1 

2 

Here,  V is  the  three-dimensional  Laplacian,  6 (x)  is  a three-dimensional 
Dirac  delta  function  and 

(1.2)  u,  (x,  X , uj)  = exp  jiii|x  - x I/cl/  (4  |x  - x |). 

i--o  ~ ~ 0 

The  signal  u^  is  incident  upon  a conve>  scatterer  (B  in  Figure  1),  giving 
rise  to  a scattered  signal  u^.  The  scattered  field  has  the  following 
integral  representation  in  terms  of  its  values  on  8 B,  the  boundary  of  B 
(See,  for  example,  BAKER  and  COPSON  [11]): 

(1.3)  Uj(x,  x^,  oO  = |ug(x',  x^,  m)  |^(^,  x',  w) 

aB 


- g (x,  x'  , 


auo 

> FrT  <!'■ 


dA. 


Here  the  integral  is  over  the  scattering  surface  in  prime  variables,  g 

is  the  Green's  function,  given  by  (1.2)  with  x replaced  by  x',  and  d/'>n 

“0 

denotes  the  outward  normal  derivative  to  aB. 
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Applying  all  of  these  results  to  (1.3)  yields 


I 


(1.6)  Uq  (x  9 X 5^) 

b - 0 - 0 


L(x  ) 


The  far  field  assumption  is  now  made,  namely  that 


(1.7)  X >>x,  x=lx|,  xonSB. 

0 


Th_.  , for  X on  3B,  |x  - x | may  be  expanded  as  follows: 


(1.8)  lx-x|=x  -x*x  + 0(x  / X ) . 

' - - 0 0 ■ 0 0 


The  result  (1.2)  is  now  inserted  in  (1.6)  and  then  from  (1.8),  a two-term 
expansion  is  used  in  the  phase  and  a one-term  expansion  is  used  in  the 


amplitude  to  yield 


(1.9)  Up  (x,  X ,u))  — 
o - -n 


exp  {2iwx  /c}  p(cox  /c) 


(4ttx  ) 

0 


Here  p((i)  x /c)  is  a phase  and  nanqe  normialiv,ed  far  field  saattering 
0 


implitude  given  by 


(1.10)  p(u)  X /c) 


/1»( 


exp  i - 2iaj  x'  • X /c}  f/dA. 


L(x  ) 


2 


I 


If  a similar  experiment  is  performed  from  the  opposite  direction, 
the  result  is 


(1.11)  p(-  CO  X /c)  = 

0 


f ||f  { ^2ia)  X'  . x^/c}}|dA. 


D(x  ) 


From  (1.10)  and  (1.11),  it  follows  that 


(1.12)  p((jj  X /c)  + p (-CO  X /c) 

0 0 


= f If  { 


2ito  x'  • X /cl  I dA 


9B 


Here  ( ) denotes  complex  conjugate. 

Applying  the  divergence  theorem  here  yields 


(1.13)  p(u)  X /c)  + p (-  O)  X /c) 
0 0 


= { 


expf-  2iut  x'  • X 


/cl  } 


dV 


= -(4iO  /c  ) y(2co  X /c)  . 

0 
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Here, 


I 


(1.14)  Y (J<)  = J exp{-  il<  • x'l  dV,  l<  = (k^,  k^,  k^). 

B 

The  function  y (K)  is  the  Fourier  transform  of  the  function 

fl,  in  B 

. 

0 X not  in  B 

This  function  is  the  characteristic  function  of  the  region  B.  If  y (x) 
is  known,  the  region  B is  known.  (Actually,  if  the  discontinuity  of 
f (x)  is  known,  B is  known.)  The  result  (1.13)  relates  the  Fourier 
transform  of  y M to  the  phase  and  range  normalized  scattering  amplitude. 
The  value  of  the  transform  variable  is 

(1.16)  k = 2w  X /c. 

~ 0 

Thus,  (1.13)  is  the  POFFIS  identity. 

It  should  be  noted  that  the  POFFIS  identity  is  a result 
derived  for  high  frequency  only.  Thus,  the  identity  is  suspect  when 
used  to  determined  y (J<)  for  small  values  of  |J<|.  Often,  as  a practical 
matter,  low  frequency  (and  hence,  small  |J<| ) information  is  simply  not 
available.  Furthermore,  the  obstacle  often  cannot  be  probed  from  all 

directions  x . Thus,  one  is  faced  with  the  limited  aperture  problem 

0 

mentioned  in  the  introduction  and  addressed  in  Section  1.3. 
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1.2  The  POFFIS  Identity  for  One-Sided  Viewinsj  of  Anomalies 
in  a Geological  Structure 


In  geophysical  or  oceanographic  exploration,  one  is  faced  with 
the  problem  of  mapping  a reflector  from  one-sided  backscatter  information, 
only.  As  an  example  of  this  type,  Figure  2 depicts  an  anticlinal 
perturbation  B of  finite  extent  in  a plane  R.  The  surface  to  be  detected 
here  is  S which  is  the  same  as  R except  on  the  anticline. 

With  Ug  denoting  the  field  scattered  by  S,  and  making  the  same 
approximations  as  in  the  previous  section,  (1.6)  is  replaced  by 


In  the  absence  of  the  anticline,  the  backscattered  field,  denoted  by  u^, 
has  the  exact  representation 


(1.18) 


Up( X 5 ^ » ^) 

K - 0 - 0 


dA. 


Thus,  the  difference  between  the  two  fields  u^  and  Up  is  given  by 


(1.19) 


■ ^R 


dA. 


Proceeding  now  as  in  the  previous  section  leads  to  the  result 

9 


'rtf. 
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(1.20) 


exp{2iw  X /c)  p(wx  /c) 
0 2 

(4Trx 

0 


with  p now  given  by 


(1.21) 


p(w  X /c)  = 
0 


J 9n 


exp{-  2iLox'  • X /c)  dA 

0 


-(4w^/c^)  Y (2wx  /c)  . 

0 


Here,  / (K)  is  again  the  Fourier  transform  given  by  (1.14,  15)  but 
with  B defined  by  Figure  2. 

In  (1.20),  Ur  is  observed  for  x having  a positive  first 

J 0 

component.  The  function  Up  does  not  have  to  be  observed  since  it 
is  a mathematical  construct  defined  by  (1.18).  Thus,  Y(ji)  is  known 
only  in  the  upper  half  k space  and  this  POFFIS  identity  does  not 

★ 

suffice  to  generate  B.  To  overcome  this  difficulty  the  region  B 

of  Figure  3 is  introduced.  This  region  is  generated  by  projecting 

★ 

B through  the  origin.  That  is,  if  ><  is  in  B,  then  -x  is  in  B . 

★ 

The  characteristic  function  of  the  region  B + B is  denoted 
by  ip(^)  ■ Thus, 


(1.22)  v(x)  = 


y(x)  X,  ^ 0 
y(-x)  X,  < 0 


Its  Fourier  transform  is  given  by 


(1.23)  .i;(k) 


/ 


B + B 


exp{ik  • x}dV  = 2 / cos(k  • x)  dV, 


B 


or,  using  (1.14) 

(1.24)  4j(Ii)  = 2 Re  Y(k). 

Therefore,  Y(]i)  ''S  determined  by  taking  twice  the  real  part  of 
y{k)  for  j<  in  the  upper  half-space  and  by  extending  this  function 
as  an  even  function  of  J<  to  obtain  its  value  in  the  lower  half- 
space. By  inverting,  one  then  finds  that 

(1.25)  y(x)  = — ^ / Re  Y(k)cos(k  • x)d^k,  Xi  > 0. 

J 

ki>  0 


In  this  case,  then  (1.20,  21,  25)  constitute  the  POFFIS  identity. 

For  a synclinal  structure  (with  B a region  below  R),  the 
result  (1.19)  is  replaced  by  its  negative.  This  change  in  sign 
persists  throughout  and  the  syncline  is  then  revealed  by  the  pro- 
cessing in  (1.25)  yielding  a function  equal  to  -1  or  zero  for 

Xi  ^ 0. 

For  a cylindrical  anomaly,  such  as  the  anticline  in 
Figure  4,  the  analysis  must  be  modified.  Firstly,  the  medium  is 
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probed  only  alone  a line  parallel  to  the  X2  axis,  say  with  X3  = 0. 
Then  in  (1.17)  and  (1.18),  X30  is  set  equal  to  zero  and  the 
integration  in  x'^  is  carried  out  by  the  method  of  stationary  phase 
[12].  The  result  is 


exp{2iwx  /c  + iir/A) 

(1.26)  Uj  - Up  = ° p(2cox  /c) 


4(2ttx  ^ 


Here, 


3/ 2 


(1.27)  p(<)  = K J exp{-iK  • X ' 1 dA;  ^ = (ki,  k2,  0) 

z 


i 


and  A is  the  cross-sectional  area  of  B in  the  (xi,  X2)  plane. 
Because  B is  a cylinder,  when  Z is  known,  B is  known. 

For  the  region  Z,  its  characteristic  function  o(xi,  X2) 
is  defined  by 


(1.28)  o(x,,  X2) 


1,  (xi , X2)  in  Z 
0 (xi , X2 ) not  in  Z 


It  then  follows  from  (1.27)  that 

(1.29)  o(k)  = p(<). 

By  using  the  same  symmetry  arguments  as  above,  one  finds  that 
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(1.30) 


a(xi,  X2) 


1 


Re  a(<)  cos(i< 


Kj  ^ 0 


2 

>^)  d K,  xi  ^ 0. 
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region  where  y(x)  = 1.  The  function  6(x-a)  has  Fourier  transform 
exp{-ika}  and  band  limited  inverse  transform 


(1.31)  6(x-a)  ~ iT"^(x-a)"^  sin[k(x-a)] 


with  k^  the  band  limits.  This  function  peaks  at  x = a,  with  peak 
value  (k^  - k_)/7T.  Thus,  the  boundaries  of  the  region  of  interest 
are  readily  recognized  from  this  band  limited  processing.  The 
purpose  of  this  section  is  to  show  that  this  is  true  in  two-  and 
three-dimensions,  as  well,  and  also  to  show  the  effects  of  aspect 
angle  limitations  in  these  cases. 

For  the  function  y{x)  defined  by  (1.15),  the  directional 
derivative  is  the  direction  defined  by  the  unit  vector  p is  given  by 


(1.32)  A(x,  p)  = -p  • Vy(x). 


If  3B  is  defined  by  = 0,  with  4>  negative  in  B and  positive 
outside  of  B,  then 

A.  /S 

(1.33)  A(^,  p)  = p • Vil>6(<l)). 


The  Fourier  transform  of  A is 


(1.34)  A(k,  p) 


= / 


p • V<J>(^)6  [<!’(£)]  exp{-ij<  • O dV. 
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Here  the  domain  of  integration  is  all  of  £-space.  The  integration 
normal  to  the  boundary  is  readily  calculated  by  exploiting  the 
properties  of  the  delta  function.  The  result  is 


(1.35)  A(k 


i.  ^)=  / 


p • n exp{-ij<  • O dA. 


Here  n is  the  outward  unit  normal  to  8B. 

The  aperture  limited  inverse  transform  of  A will  be  denoted  by 
Ai.  If  Ap(l<)  denotes  the  aperture  in  space,  then 


(1.36)  Ai(x,  p)  = 


dA  p • n expiik  • (x  - ') 


Ap(k)  3B 


The  apertures  of  interest  are  those  for  which  k = |J<|  ranges  between 
two  values  k and  k^  and  the  angle  of  k is  restricted  to  some  region  fi. 


Thus , 


(1.37)  A,(x,  p)  = 


k^. 


dA  p • n exp{ i k • (x  - ^) 


SI  9B 


Here  dfi  is  the  solid  angle  element  in  k space. 

The  four  fold  integral  over  SI  and  3B  can  be  calculated  by  the 
method  of  four-dimensional  stationary  phase  [12].  The  integration  over 
k can  then  be  done  explicitly.  The  method  is  described  in  [10],  but 
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the  form  of  the  result  presented  here  is  due  to  Armstrong  [14].  The 
conditions  that  the  phase  be  stationary  are  as  follows.  Firstly,  the 
point  ^ on  3B  must  be  such  that  - £ is  perpendicular  to  9B.  Secondly, 
the  angles  of  k are  such  that  ^ - K and  l<  are  colinear  or  anticol inear. 
Thus,  the  stationary  points  in  ^ and  J<  are  such  that  ^ ^ and  k lie 

along  the  normal  to  9B  at  the  stationary  point. 

For  each  choice  of  x there  may  be  one  or  more  stationary 
points  (or  none  at  all!).  If  the  contribution  from  each  stationary 
point  is  denoted  by  A2(x_,  p),  then  one  finds  that  asymptotically 

(1.38)  A2(x.  P)  = (2it)'^  n*p  (1  - mKiD)'''^(1  - yK2D)"^^ 

• I sin(kD)  - i[cos(kD)  ’ 1]  | 

k = k_ 

Here, 

(1.39)  D = |x  - i|  ; 

Ki  and  K2  are  the  principal  curvatures  of  3B  at  the  stationary 
point;  u = +1  or  -1  according  to  whether  £ - ^ and  n are  anticolinear 
or  colinear  at  the  stationary  point. 

The  real  part  A2  is  seen  to  have  the  same  qualitative 
behavior  for  D small  as  the  one-dimensional  result  (1.31). 

Furthermore, 
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(1.40)  lim  A2(x>  p)  = irp  (k  - k ). 

D - 0 

For  each  point  on  9B  the  contributions  from  stationary  points 
with  D = 0 will  dominate  all  other  contributions  to  Ai(>^,  p)  if 
the  bandwidth  is  sufficiently  broad.  There  can  be  at  most  two 
such  cont 'ibutions.  This  will  occur  if  the  region  contains  both 
the  vectors  k = n and  k = -n  for  the  given  stationary  point, 

In  [9],  the  following  test  of  this  method,  as  applied  to 
the  POFFIS  identity,  was  carried  out.  The  exact  backscatter  solu- 
tion for  a sphere  of  unit  radius  was  used  to  generate  p(aiXo/c)  for 
the  POFFIS  identity  (1.13).  The  angular  aperture  was  taken  to 
be  unrestricted;  i.e.,  all  angles  were  used.  The  vector  p = (1,  0,  0) 
was  used  and  Ai(x,  p)  was  calculated  for  x?.  = X3  = 0.  The  angular 
integrations  were  carried  out  numerically  and  the  integration  in 
k was  done  numerically  with  the  trapezoid  rule.  An  example  of 
the  output  of  this  processing  is  shown  in  Figure  5.  The  circles 
represent  the  theoretical  results  based  on  the  asymptotic  analysis 
presented  here.  In  this  case,  k_  = 15.75  : 5tt,  k^  = 31.5  ~ IOti. 

At  the  peak,  the  percentage  error  between  theory  and  computation 
was  .3%. 

The  two-dimensional  analog  of  the  results  presented  here 
is  carried  out  in  exactly  the  same  manner.  In  (1.37),  9B  would 
now  be  a curve  and  dfi  is  replaced  by  do,  the  polar  angle  of 
= (ki,  k?).  The  factor  (2tt)"’  is  replaced  by  (Ptt)'^.  The  result 
(1.38)  is  now  replaced  by 
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(1.41) 


A (x^,  p)  = (2it)”*  n*p  (1  - uKD)"'^ 

• D‘*  I sin(kD)  -i[cos(kD)  - 1]  ^ 


Here  y and  D are  as  defined  below  (1.38). 

In  [15],  the  two-dimensional  result  was  tested  for  an 
ellipse  with  semi -axes  of  length  one  and  two.  The  result  for  full 
angular  aperture  k range  9 to  27  and  p = (1,  0)  is  shown  in 
Figure  6.  The  third  dimension  (the  value  of  A2)  is  laid  back  down 
in  the  plane.  On  each  line  the  height  is  normalized  with  respect 
to  maximum  height.  Furthermore,  the  maximum  on  each  line  is  tested 
against  the  absolute  maximum.  When  this  relative  maximum  falls  below 
a critical  value,  the  entire  line  is  zeroed  out.  This  insures, 
for  example,  that  "noise"  will  not  be  enhanced  outside  of  the 
vertical  extent  of  the  ellipse.  Nonetheless,  near  the  top  and 

/V 

bottom  of  the  ellipse,  where  p . n is  nearly  zero,  there  is  clearly 
some  loss  in  resolution.  This  region  could  be  resolved  better  if 
p were  chosen  to  be  (0,  1). 

In  Figure  7,  limited  aperture  process  is  shown  for  Q 
consisting  of  two  quadrants.  There  is  some  "spill-over"  outside  of 
the  angular  aperture,  but,  for  the  most  part  the  ellipse  is 
reproduced  as  in  the  previous  diagram,  but  only  for  the  region  for 
which  the  normals  to  the  ellipse  lie  in  the  angular  aperture  in 


J<  space. 

1 
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In  [10],  further  examples  are  carried  out  in  which 
numerical  data  is  generated  for  backscattering  by  a circular 
cylinder.  The  results  confirm  the  POFFIS  method  of  Section  1 as 
well  as  the  Fourier  analysis  of  this  section. 

It  should  be  noted  that  near  the  evolute  of  3B,  (1.38) 
or  (1.41)  becomes  invalid.  In  [14],  ARMSTRONG  shows  that  the  con- 
tribution from  that  region  is  extremely  small  relative  to  the  peak 
value  of  A?  at  D = 0. 


i 
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1.4  POFFIS  in  the  Time  Domain 


Returning  to  (1.9),  the  inverse  time  transforms  of  and  P 
are  defined  to  be  and  V,  respectively.  From  (1.9)  it  follows 
that 


(1.42)  V(xo,  t)  = (4ttxo)'  U.(x  , t + 2xo/c). 

o - 0 


That  is,  the  Fourier  transform  of  the  phase  and  range  normalized 
scattering  amplitude  is  a range  normalized  and  time-delayed 
scattering  amplitude  in  the  time  domain  and  thus,  is  as  easily 
observable  as  p(wxo/c).  The  objective  of  this  section  is  to  show 
how  to  use  this  result  and  the  basic  POFFIS  identity  to  derive  a 
POFFIS  identity  in  the  time  domain. 

The  Fourier  inversion  in  (1.13)  is  first  expressed  as  an 
integral  in  polar  coordinates.  The  result  is 


(1.43)  y(x)  = 


f f 

4tt3c  J J 


dfi{p((JXo/c)  + P (-wxo/c)} 


• exp{2icoXjj  • >^/c} 


Here,  12  denotes  the  unit  sphere  with  variable  xo  and  differential 

solid  angle  element  on  it,  dJ2.  The  volume  element  in  polar 

2 

coordinates  has  a factor  of  w ; however,  the  integrand  also  has  a 
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(1.46)  y(x)  = - 


00 

— / duj  / dn| 

l^'cj  J ' 


p(u)x^/c)  exp{2i(i)x^*  x/cl 


/\  \ 
+ p(-(jjx^/c)  exp{-2i(iJx^  • x/c}  / 


Here,  p in  the  first  line  is  defined  by  (1.10)  and  in  the  second 
line  by  (1.11).  However,  again  using  the  fact  that  D(xJ  = L(-x^), 
the  result  can  again  be  written  as  an  integral  over  a sphere: 


(1.47)  y(x)  = - 


00 

I 

-0°  n 


dfi  p(<x)x  /c)  exp{2iwx  • x/c} 
0 0 - 


The  inverse  time  transform  of  p was  defined  to  be  V. 
Therefore,  (1.47)  can  be  rewritten  as 


(1.48)  y(x)  = - 


kf 


d^l  \/(x  , -2  X • x/c) 


This  is  the  POFFIS  identity  in  the  time  domain. 

The  result  (1.47)  also  provides  an  alternative  to  the 
POFFIS  identity  (1.13).  In  the  form  (1.46),  the  requirement  of 
combining  "front  side"  and  "back  side"  observations  before  Fourier 
synthesis  has  been  dispensed  with.  There  is  still  a double  covering 
in  space  in  (1.47),  but  this  may  be  overcome  by  using  (1.44)  to 
conclude  that 
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(1.49)  ,(x)  = 


Re 

2n^c 


i-f 


dQ  p(ujX  /c)expi2i..x  • x/c}  . 


In  this  form,  the  integral  is  again  a Fourier  transform.  When  (1.9) 
is  used  and  the  integral  is  expressed  again  in  terms  of 


(1.50) 


k = (.X  /c  , k = X , 

“0  n 


the  result  becomes 


(1.51)  ,(x)  = - 8ti’ ' Re  / k’  d ’k  Xj,u^.(kx  , kx  , ck) 


• expi2ik*x  - 2ikx  I 

— — Q 


Similarly,  one  can  show  that  for  j = 1,  2,  3, 


:)y(x)  f _ ^ 

(1.52)  — = 16ir  ' Im  / k-k  ^d  ^k  x.  U(.(kx  , kx  , ck) 

3Xj  J ^ 0 b 0 


• exp{?ik*x  - 2ikx 


In  these  two  equations  x,  may  be  a function  of  x^  (i.e.,  the  ob- 
servation distance  may  be  different  in  different  directions). 
Thus,  in  general,  x must  remain  under  the  integral  sign. 
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2. 


An  Inverse  Method  for  Determining  Small  Inhomoqenei ties 
In  a Medium 


When  a medium  is  known  to  have  small  and  perhaps  slowly 
varying  inhomogeneities,  the  method  of  Section  1 is  inapplicable  to 
the  problem  of  reconstructing  the  medium.  The  method  presented 
here  is  a procedure  for  determining  such  variations.  This  method 
would  be  appropriate  to  the  detection  of  velocity  gradients  in  the 
seabed  as  well  as  to  the  general  problems  of  subsurface  mapping  or 
geophysical  exploration.  Mathematically,  the  inverse  problem 
considered  here  is  to  determine  coefficients  in  a system  of 
equations  governing  a wave  propagation  problem.  It  is  assumed  that 
the  reference  values  of  these  coefficients  (often  constants)  are 
available  and  that  variations  from  these  reference  values  are  to  be 
determined  from  observations  of  the  field  arising  from  known  input 
sources.  Physically,  the  coefficients  to  be  determined,  usually 
characterize  the  medium  velocities  or  the  acoustic  impedance  or 
some  similar  quantity. 

For  clarity  of  exposition,  the  discussion  here  will  be 
restricted  to  media  in  which  propagation  is  governed  by  the  three- 
dimensional  wave  equation.  In  reference  [16]  by  the  authors,  the 
inverse  problem  for  electromagnetic  waves  and  elastic  waves  is 
discussed,  as  wel 1 . 

An  essential  feature  of  the  inverse  method  presented  here 
is  that  an  integral  equation  is  derived  for  a function  which 
characterizes  the  velocity  or  impedance  variation.  This  equation 


i 


is  a Fredholm  integral  equation  of  the  first  kind;  it  has  parameters 
in  it  which  characterize  the  source  and  receiver  locations.  In  a 
number  of  cases  of  source-receiver  configurations  of  practical 
interest,  this  integral  equation  is  invertible.  Thus,  a velocity  or 
impedance  profile  (depth  section)  is  obtained  by  direct  processing 
of  the  observed  data  (time  section)  itself;  i.e.,  by  performing 
weighted  quadratures  on  the  data.  This  direct  inversion  of  time 
section  to  depth  section  involves  only  one  theoretical  assumption: 
The  subsurface  variation  must  be  "small."  Even  this  limitation  is 
not  overly  restrictive  as  can  be  seen  from  one  of  the  examples  below 
in  which  a 20X  velocity  variation  was  successfully  migrated.  In 
fact,  the  real  world  data  restrictions--noise , attenuation,  dis- 
cretization and  finiteness  of  observations,  etc. , --are  usually  of 
greater  concern  than  this  theoretical  linearization  assumption. 
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2.1  An  Integral  Equation  for  Three-Dimensional 
Velocity  Variation 

It  is  assumed  here  that  three-dimensional  wave  propagation 
is  governed  by  the  equation 

(2.1)  [V^  - v'Mx,  z)9^^]  U(x,  z,  t)  = 0;  x = (x,y). 

Waves  propagate  in  the  semi-infinite  medium,  z 0,  measured  down- 
ward. The  operator  9^  denotes  a partial  derivative  with  respect 
to  time.  It  is  the  coefficient  v(>^,  z)  which  is  to  be  determined 
from  probes  and  observations  made  at  the  upper  surface,  z = 0.  For 
many  cases  of  interest,  the  surface  z = 0 separates  two  media  of 
greatly  differing  impedances  (e.g.,  air-water,  air-earth).  Hence, 
an  appropriate  boundary  condition  for  (2.1)  which  introduces  the 
probes  is 

(2.2)  = 't>(t,  x;  A),  z = 0,  A = (^H)  . 

The  sources  ‘I'  characterized  by  the  vector  ^ may  be  many  types  (See 
[16]  for  a discussion  of  plane  wave  sources).  For  definiteness,  here, 
the  discussion  will  be  restricted  to  impulsive  point  sources  located  at 
>^  = _A,  z = 0.  Thus,  the  boundary  condition  actually  treated  is 

(2.3)  9^U  = -6(t)  6(x  - A),  z = 0 . 
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To  complete  the  specification  of  the  direct  problem  for  U,  it  is 
assumed  that 

(2.4)  U 0,  t < 0. 

The  reference  value  of  the  velocity,  which  is  assumed  known,  will 
be  denoted  by  c(^,  z).  Then,  the  variation,  denoted  by  i(x,  z) 
is  defined  by  the  equation 

(2.5)  v'^(x,  z)  = c'^(x,  z)[  1 + .(x,  z)]  . 

It  is  the  variation  i which  is  to  be  determined  by  means  of 

observations  of  U made  at  surface  points  x = j_,  z = 0. 

The  total  field  U is  decomposed  into  a primary  field  Uj, 
which  is  the  solution  in  the  absence  of  the  variation  «,  and  a 
scattered  field  U^,  which  is  the  response  to  a.  Thus 

(2.6)  U(t,  X,  z;  a)  = Uj(t,  X,  z;  A)  + Ll^(t,  x,  z;  ^) , 
where  Uj  satisfies  the  following: 

(2.7)  [V  - ] Uj  = 0,  z 0; 

(2.8)  Uj  = --(t)  (x  - ■),  z = 0; 

(2.9)  Uj  0,  t 0. 
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The  argument  ^ is  introduced  into  U,  Uj,  and  to  emphasize 
their  dependence  on  the  source  location.  It  then  follows  from 
the  problems  for  L)  and  Uj  that  satisfies  the  following: 

(2.10)  [v^  - U3  = -a  c'^3^^  U; 

(2.11)  9^03  =0,  z = 0; 

(2.12)  U3  = 0,  t - 0. 

An  adjoint  function  V is  now  introduced,  satisfying  the 

following: 

(2.13)  (V^  - c""9^")  V = 0,  z > 0; 

(2.14)  9^V  = -(T  - t)  H(t  - t)  6(x  - O,  z = 0;  i = (S,n); 

(2.15)  V = 0,  t > T. 

It  will  be  seen  below  that  C,  t denote  the  locations  and  times  of 
observations  at  the  upper  surface.  Furthermore,  by  comparing  the 
problems  for  Uj  and  V,  it  can  be  seen  that 

(2.16)  9^'  V(t,  X,  z;  C)  = U(t  - t,  i,  z;  x ) . 
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Upon  applying  Green's  theorem  (in  space  and  time)  to  the 


quantity. 


the  following  integral  equation  is  obtained: 


(2.17) 


rxj  oo 

/dz  /c 


I 

:x,  z)Jd 


dxdy  «(x,  z)  c'  (x,  z)/dt  U(t,  x,  z;  A)U,(''-t,  , z;  x) 


dtU3(t,  f,  0;  A)( I - t). 


The  right  side  here  is  a function  of  the  field  observations  at 
the  upper  surface  and,  hence,  known.  The  left  side  has  two  unknowns, 
namely,  and  U^.  (A  nonlinear  system  for  these  two  unknowns  is 
supplied  by  (2.17)  and  (2.10,  2.12).)  However,  appears  only  in 
U and,  therefore,  only  through  the  product  . However,  from 
(2.10)  it  is  seen  that  is  itself  of  the  order  of  . Thus, 
for  ' small,  it  is  expected  that  O'  U can  be  reasonably 
approximated  by  Uj.  In  this  case,  (2.17)  becomes  an  integral 
equation  for  ‘ alone,  namely. 


( 


00  OO  1 

(2.18)  /dz  / dxdy  a(x,  z)c"^()^,  z) 

0 -OO  0 

= j dt  U3(t,  i,  0,  A)(t-  t) 
0 

This  is  the  basic  linear  integral  equation  for  a{x,  z 


Jdt  Uj (t,  X, 
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z;  A)Uj(T-t,  z;  x) 


). 


2.2  Direct  Inversion  for  Backscatter  Over  a Medium  with 
Two-Dimensional  Velocity  Variation 

The  integral  equation  (2.18)  will  now  be  specialized  as 
follows.  The  reference  velocity  c will  be  taken  to  be  a constant, 
so  that  Uj,  defined  by  (2.7  - 9)  is  given  by 

(2.19)  Uj(t,  X,  z;  A)  = 6(t  - |x  - A|/c)/(2tt|x  - 2^|  ). 

The  variation  in  a will  be  assumed  to  be  independent  of  y,  so  that 
a(x,  z)  will  be  replaced  by  a(x,  z).  In  this  case,  both  the  sources 
and  receivers  will  be  restricted  to  the  x axis.  Finally,  only 
backscattered  (or  CDP  stacked)  observations  will  be  made,  so  that 

(2.20)  A = i = (T,  0). 


For  this  case,  it  is  possible  to  carry  out  the  y and  t 
integrations  in  (2.18)  to  obtain 

oo 

(2.21)  j dz 
0 


/ 


dx  u(x,z)H(ct/2  - p)  [(ct/2)^-  p^] 


= -21(710) 


I 


dt 


Us(t, 


C,  z;  ^)(t-  t). 


Here, 


1 
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(2.22)  p ■ [ (x  - {)“  * 2’] 


This  integral  equation  can  be  solved  by  transform  techniques. 
The  result  is 


(2.23)  a(x,  z)  = 2ic^7T 


OO  OO  CO  CO 


k (I'-Tt) 


“OO  — OO  *“0O 


0 0 


• U<;(t,  C.  z;  0 exp{2ik  (x-C)  - 2ik  z + iu)T}>; 

J “ “ 1 3 ) 

0)  = c[sgnk  ] (k  ^ + k 

3 1 3 

This  is  the  direct  inversion  formula  for  the  prescribed  source 
receiver  configuration  for  this  section. 

In  [18],  by  the  authors,  this  result  was  implemented  for 
high  frequency  data  synthetically  produced  for  a number  of  layered 
models.  Because  of  the  high  frequency  bandlimiting,  the  data  was 
processed  for  9oi  /9z  = a'  rather  than  for  a(x,  z),  itself.  This 
merely  requires  multiplication  of  the  integrand  in  (2.23)  by  -2ik^. 
In  accordance  with  the  results  of  Section  1.3,  not  only  the 
location  of  discontinuities  may  be  determined  but  the  magnitude 
of  the  discontinuities  may  be  determined,  as  well.  This  direct 
inversion  procedure  produces  a mapping  of  the  interfaces  and  an 
estimate  of  the  velocity  variations  across  them. 
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The  model  used  embodies  other  real  world  restrictions  band 
limiting,  namely; 

(i)  the  observations  are  made  only  at  discrete 
points  on  the  line; 

(ii)  the  observations  are  made  only  over  a 
finite  length  of  the  line. 

In  these  examples,  the  spacing  between  r points  is  100  feet,  the 
maximum  array  length  is  8,000  feet  and  the  bandwidth  is  6 to  24 
Hertz. 

Synthetic  data  was  generated  for  tilted  planes  at  angles 
up  to  75°.  Direct  inversion  reproduced  these  planes  with  three- 
place  accuracy  both  in  tilt  angle  and  velocity  variation. 

In  an  earlier  paper  by  the  authors  [19],  a direct 
inversion  formula  was  presented  for  the  parabolic  approximation  to 
the  wave  equation  used  by  CLAERBOUT  [17]  for  15°  wave  equation 
migration.  Synthetic  data  for  tilted  planes  is  used  and  the 
direct  inversion  is  carried  out  analytically.  The  results  show  an 
error  in  tilt  angle  and  in  velocity  estimation  on  the  order  of  the 
fourth  power  of  the  angle  of  inclination.  Numerical  tests  confirm 
this  result,  as  well.  This  kind  of  error  is  already  known  to 
users  of  15°  wave  equation  migration.  The  error  is  due  to  the 
underlying  parabolic  approximation  and  not  a flaw  of  the  direct  in- 
version procedure  or  the  wave  equation  migration  procedure.  Since 
the  computational  requirements  of  direct  inversion  for  the  wave 
equation  are  no  more  difficult  than  those  for  the  parabolic 
approximation,  the  former  is  preferred  to  the  latter  by  the  authors. 
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Figure  8 shows  a subsurface  configuration  for  which 
synthetic  data  was  graciously  provided  to  us  by  the  research  group 
at  Marathon  Oil.  The  time  section  provided  by  them  is  depicted  in 
Figure  9;  Figure  10  is  the  result  of  our  direct  inversion  procedure 
and  Figure  11  shows  our  estimate  from  the  output  of  various 
relevant  quantities  in  the  model.  The  lower  two  sets  of  numbers 
exhibit  errors  of  less  than  4%,  while  above  the  level,  the  errors 
are  less  than  1%. 

The  program  was  implemented  on  a time-sharing  system  on 
a Burroughs  6700  at  the  University  of  Denver.  Output  data  was 
generated  at  a rate  of  5 mi li seconds  per  record.  The  number  of 
records  is  defined  as  the  product  of 

(No.  of  geophones)  x (No.  of  subsurface  output 
points)  X (Average  number  of  non-zero  time 
samples/trace) . 

Modifications  of  the  original  program  have  now  improved 
that  rate  to  2.4  mill  seconds  per  record. 
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OO  00 


00  OO  J 


(2.26)  a(x,  z)  = 2ic^7T'" 


J d^dn  y* dk  dk^y  dk  JctiJcltj 


k (i^-it) 


— OO  -OO 


0 0 0 


• Uc.(t,  C.  z;  C)  exp{2i[k  (x-C)  + k (y-n)-k  z]+v^T}i; 

^ ~ — 1 2 3 f 


oj  = c[sgnk  ] (k^  + k^  + k^)^ 

3 12  3 


A computer  program  could  be  developed  to  implement  this 
formula  much  as  was  described  in  the  previous  section.  Output  data 
would  be  generated  at  the  same  rate  per  record  discussed  at  the 
end  of  the  section.  However,  the  number  of  records  increases  by  a 
factor  (say,  N)  for  the  new  dimension  in  the  geophone  array  and  by 
a factor  (say,  N)  for  the  output  array.  Thus,  the  computer  time 
would  increase  by  a factor  of  N^.  Thus,  for  the  present,  only 
analytical  verifications  have  been  carried  out. 
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2.4 


Direct  Inversion  for  a Case  with  Separated  Source 
and  Receiver 


Here  it  is  assumed  that  = c<(z),  i.e.,  thata  varies  with 
z alone.  Intuitively  then,  one  would  expect  that  only  one  experiment 
would  be  necessary  to  determine  a.  Again,  it  is  assumed  that  c is 
constant  so  that  Uj  is  given  by  (2.19).  In  this  case,  ^ and  £ are 
taken  to  be  fixed  with 


(2.27)  C = - i . 


In  (2.18),  the  x,  y and  t integratiotis  can  all  be  performed 
yielding  an  elementary  integral  equation  which  can  readily  be 
inverted  to  yield 


2(A^  + z")7c 


(2.28)  a(z)  = -4t,c  dt  \ ~ 


f 


2 + 2z^ 


0 


^ (A^+z^)'^ 


U5(t,  A,  0;  -2) 


For  backscatter,  it  is  only  necessary  to  set  A = 0 here  to  obtain 
the  result 

2z/c 

(2.29)  ‘(z)  = -4nc  j dt  [4z/c  - t]  Ug(t,  0,  0;  0). 
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2.5  Direct  Inversion  for  a One-Dimensional  Problem 


In  one  dimension,  one  can  generate  synthetic  wide  band 
data  by  straightforward  and  economic  means.  This  was  done  by 
GRAY  [20]  as  a prelude  to  the  analysis  of  strongly  depth  dependent 
velocities  in  three  dimensions.  An  example  of  the  output  of  this 
analysis  is  shown  in  Figure  12.  The  solid  line  is  the  assumed 
a (z)  while  the  dots  are  the  result  of  direct  inversion  on 
synthetically  generated  wide  band  backscattered  data. 

In  [16],  the  authors  also  treated  a case  in  which  c(x) 
was  not  constant.  The  function  Uj,  the  time  transform  of  Uj,  was 
expressible  as  a sum  of  exponentials  and  the  integral  equation  was 
still  invertible  in  closed  form. 
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2.6  Direct  Inversion  In  Free  Space 


Often  in  oceanographic  and  seismic  exploration,  the  source 

is  placed  below  the  surface  z = 0.  In  this  case,  the  basic  problem 

(2.1,  3,4)  should  be  replaced  by  a source  problem,  with  the  delta 

function  moved  from  the  right  side  of  (2.3)  to  the  right  side  of 

(2.1)  and  the  boundary  condition  (2.3)  imposed  on  a new  interface 

above  z = 0,  say  at  z = z <0. 

0 

If  the  effects  of  the  reflections  from  the  interface  at 

z can  be  accounted  for  (by  no  means  an  easy  task),  then  such  an 
0 

experiment  can  be  modeled  by  a free  space  problem  in  which  the 
medium  is  "known"  for  z < 0 and  assumed  to  vary  only  for  z > 0. 

In  this  case,  one  still  obtains  (2.18)  as  the  basic  integral 
equation,  but  with  Uj  now  a solution  of  the  corresponding  source 
problem  in  free  space.  For  c = constant,  the  effect  of  this  modifi- 
cation on  Uj  is  to  replace  the  factor  of  2 on  the  right  side  of 
(2.19)  by  a 4.  Since  two  factors  of  Uj  appear  in  the  kernel  of  the 
integral  equation  (2.18),  the  results  (2.23,  26,  28)  need  only  be 
modified  by  the  introduction  of  a multiplier  of  4 on  the  right 
side. 


T 
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3.  Non-Uniqueness  in  the  Inverse  Source  Problem 

In  the  previous  sections  the  presentation  assumes  that 
the  probing  signal  is  an  impulsive  point  source  (See  1.1,  2.3).  In 
practice,  this  is  never  the  case.  Thus,  one  is  faced  with  the  prob- 
lem of  "stripping  away"  source  effects,  i.e.,  performing  source 
deconvolution,  before  proceeding  to  the  inversion  problem  at  hand. 

In  other  problems,  the  source  itself  is  the  ultimate  goal 
of  the  inversion  process.  This  would  be  the  case  either  for  detec- 
tion problems  or  for  source  synthesis  problems  (such  as  antenna 
design  problems)  in  which  the  objective  is  to  create  a source  to 
produce  a given  radiated  field. 

The  objective  of  this  section  is  to  address  the  inverse 
source  problem.  The  presentation  follows  the  lines  of  [21]  by  the 
authors.  It  will  be  shown  that  there  is  a great  deal  of  non- 
uniqueness in  this  problem.  That  is,  a source  cannot  be  completely 
reconstructed  from  observations  of  the  field  it  radiates.  The 
reason  is  that,  in  general,  a certain  part  of  the  source  simply 
does  not  radiate  and,  hence  yields  no  clue,  in  the  radiated  field, 
to  its  presence.  Thus,  the  solution  to  the  inverse  source  problem  is 
non-unique.  Analytic  characterizations  of  this  non-uniqueness  will 
be  stated. 

For  source  deconvolution,  the  non-uniqueness  is  not  as 
serious  a problem  as  it  might  at  first  appear  to  be.  The  reason  is, 
j that  only  the  radiating  part  of  the  source  is  of  interest  in  these 


i 

. « ‘ ^ . f.  . 


1 

problems.  This  is  so  because  only  the  radiating  part  of  the  source 

can  effect  the  scattering  obstacle  or  medium.  I 

i 
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3.1  Analysis  of  the  Direct  Problem  for  Acoustic  Waves 

It  is  assumed  that  U(x,  t)  is  a solution  of  the  following 
problem  in  free  space: 


(3.1)  [v^  - c’'9^2]  U(x,  t)  = -F(x,  t); 


(3.2)  UHO,  F=0,  t<t. 

0 


The  source  F is  assumed  to  be  confined  to  some  sphere  x a,  to  be 
denoted  by  D.  This  region  is  itself  confined  to  the  interior  of 
some  larger  region  Dg.  Observations  of  the  radiated  field  will  be 
made  on  90,  the  boundary  of  D. 

The  Fourier  time  transform  and  space-time  transforms 
are  defined  by  the  following: 


(3.3) 


a 

u(x,  U))  = jd 


dt  U(x,  t)  exp{iwt}  ; 


(3.4)  u(j<,  oj) 


00  oo 

hi 


•00  -00 


dV  U(x,  t)  exp{-ij^  • x + iwt} 


00 

J dV  u(x,  to)  exp{-ik  • x}  . 

— 00 
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The  time  reduced  problem  equivalent  to  (3.1,  2)  is 


(3.5)  (v^  +0)2/0)  u(>L.  to)  = -f(x,  to) 

with  u outgoing, 

(3.6)  u(x,  0))'^-  u (x,  to)  exp{  ito  x/c}  /(4tt  x),  x . 

~ 0 

The  solution  to  this  problem  can  be  expressed  as  a con- 
volution of  the  source  with  the  free  space  Green's  function: 

(3.7)  u(£,  0))  = / dV  f (x ' , to)  g(r,  oi) . 

X < a 

Here,  dV  is  the  differential  volume  element  in  the  prime  variables, 

(3.8)  r = |x  - X ' i . 

and 

(3.9)  g(r,  to)  = exp{iwr/c}  /(4nr). 

The  solution  on  90  is  of  interest.  The  origin  of  the 
coordinate  system  is  taken  to  be  in  D.  Then  for  x on  90,  x > a. 

In  this  case,  using  the  spherical  harmonic  expansion  for  g 
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(JACKSON,  [22],  p.  742),  the  solution  u may  be  expressed  as  follows: 


(3.10) 


u(x,  u) 


00  m 

iwc”' XI  £ 

m=0  n=  -m 


mn 


Here, 


(3.11) 


^mn 


a 

j x'^dx’, 

0 


j and  h are  the  spherical  Bessel  function  and  Hankel  function 
'^m  m 

of  the  first  kind,  respectively,  and  is  the  spherical  harmonic 
of  order  mn.  The  functions  f^^  are  defined  by 

IT  ZtT 


(3.12) 


(d) 


j sine  de  J 


d(!)  f(x,  w)Y|^^  (e,<))) 


These  functions  are  the  coefficients  of  f in  its  spherical  harmonic 
expansion 


(3.13)  f(x. 


w)  = 1] 

m=0 


m 

E 

n=  -m 


mn  mn 


Since  the  spherical  harmonics  are  a complete  set  of 
functions,  knowledge  of  the  coefficients  f^^  constitutes  knowledge 
of  f itself.  However,  from  (3.10,  11)  it  is  seen  that  the  radiated 
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field  (i.e.,  u for  x > a)  is  a function,  not  of  the  but  only 

of  their  projections  on  the  Bessel  functions  Stated  another 
way,  the  radiated  field  depends  only  on  the  projections  of  the 
source  on  the  "doubly"  infinite  set  of  functions  {j  (wr/c)Y  (0,(ti)}. 
Since  this  class  of  functions  is  known  not  to  be  complete,  knowledge 
of  these  projections  does  not  suffice  to  determine  f uniquely. 

To  demonstrate  this  non-uniqueness,  one  need  only  produce 
a non-zero  function  f for  which  all  of  the  projections  c in  (3.11) 
are  zero.  Such  an  example  is  provided  by  the  function 


In  this  example,  all  of  the  functions  f , m > 0,  are  zero,  while 
f is  not.  However,  the  projection  of  f on  j is  zero.  Thus 

00  y r 00  0 

for  (3.10),  the  radiated  field  is  zero,  while  the  source  is  not. 
The  inverse  time  transform  in  (3.14)  produces  a source  in  space- 
time  confined  spatially  to  the  region  r < a and  temporally  to  the 
interval  -a/c  < t < a/c. 

An  alternative  characterization  of  this  non-uniqueness 
can  be  deduced  as  follows.  One  solves  the  space  time  transformed 
problem  for  u and  inverts  the  time  transform  to  obtain  the 
representation  of  the  solution 


! 

I 

i 
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00 


(3.15)  U(£,  t)  =Tr'*c/16  2]+  J k"'d^k  f(J<,  + ck)exp{ij<  • x + ckt}  . 

-00 

Thus,  the  solution  is  seen  to  depend  upon  f(j<,  co)  only  on  the  two 
sheets  of  a hyper-cone  where  k = lw|  /c.  Thus,  the  source  need  only 
have  Fourier  transform  equal  to  zero  on  this  hyper-cone  in  order 
that  the  radiated  field  be  zero. 

As  an  example,  the  source  (3.14)  has  Fourier  transform 
a /a 

(3.16)  f(J<,  w)  = 1 - J"  j^(wx/c)j^(kx/c)  J ^ {^x/c)x^dx 

0 ' 0 

which  indeed  vanishes  when  w = +ck,  but  is  clearly  not  identically 
zero. 

For  the  electromagnetic  case  with  current  source  distribution 
J()^,  t),  there  will  be  no  radiated  field  if 

(3.17)  [I  - kk]  i(k,  + ck)  = 0. 

Here  I is  the  identity  operator  and  kk  is  the  dyadic  operator  that 
yields  the  radial  part  of  j in  J^-space.  The  operator  in  (3.17) 
therefore  produces  the  traneveree  ^as  opposed  to  radial)  part  of  j. 
Thus,  any  j which  is  purely  radial  will  produce  a non-radiating  field. 
In  addition,  there  will  be  no  radiated  field  if  only  the  transverse 
part  of  j vanishes  on  the  appropriate  hyper-cone. 


1 
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3.2  Far  Field  Analysis  of  the  Direct  Problem 

In  the  far  field,  x > > x',  the  Green's  function  in  (3.7) 
can  be  expanded  just  as  in  Section  1 (see  equation  1.8).  Then  one 
finds 

(3.18)  u (x,  w)  = / dV  f(x,  a))exp{  -itox  • x/c} 

“ J 

X £ a 

= f(a)X/c,  0)). 

Here  u is  the  phase  and  range  normalized  far  field  scattering 
0 

amplitude  defined  by  (3.6).  Thus,  observation  of  the  far  field  at 
all  frequencies  and  in  all  directions  provides  an  asymptotic 
expression  for  that  part  of  f(l<,  w)  which  produces  the  radiated 
field. 

It  should  be  noted  here  that  (3.18)  provides  an  asymptotic 
solution  of  the  inverse  source  problem  for  f on  the  appropriate 
hyper-cone. 
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3 . 3 The  Inverse  Source  Problem 


The  objective  here  is  to  develop  methods  for  obtaining 
information  about  a source  distribution  solely  in  terms  of  obser- 
vations made  on  9D. 

The  function  v(^,  w)  is  introduced,  denoting  any  solution 
of  the  homogeneous  reduced  wave  equation  in  the  domain  D.  Then 
Green's  theorem  is  applied  on  the  domain  D to  the  quantity 
v(>^,  co)  [V^  + 0)^  c u()^,  (jj)  - u(>^,  co)  [V^  + co^c"^]  \/{x,  (jj) . 

Here,  u is  a solution  of  (3.5).  The  result  is 

(3.19)  J v(x,  w)  f(x,  u))  dV  = y|u  - V 1^1  dA. 

X < a 3D 

For  any  choice  of  v(£,  co)  of  the  prescribed  type,  this  is 
ar  ’-'cegral  equation  for  f()^,  w).  A particular  choice  suggested  by 
the  disc.us:ion  of  the  previous  sections  is  any  function  in  the  class 

(3.20)  = j_„(a)x/c)  Y (e,())),  m = 1,  2 |n|  < m 

' ' mn  ^mn  mn  ' i i _ 

with  V = (3.19)  becomes, 

it 

iu  _EL  . 4/*  Ml  dA. 

» 3n  3n^ 

3D 

Here,  the  c„„'s  are  defined  by  (3.11).  The  source 
mn 
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(3.22) 


f (>i,  (Jj)  = 
1 


mn 


mn 


X < a 


; X > a 


4- 

will  produce  the  same  radiated  field  as  f itself.'  Thus,  f (x,  w) 

1 ~ 

would  serve  as  an  equivalent  source  for  the  purpose  of  source 
deconvolution. 

Another  choice  of  v(>^,  w)  is  the  plane  wave 

(3.23)  v(x,  0),  k)  = exp{i  j<  • >^  } , k = |aj|/c. 

Here  k may  range  over  any  of  the  continuum  of  directions  in  three 
space.  Now  (3.19)  becomes 


(3.24)  f(k,  0)) 


^ ij£  n • u - 


8n 


exp{ij<  • 21}  dA,  k = |w|/c. 


8D 


This  is  the  exact  result  for  which  (3.18)  is  the  far  field  asymptotic 
expansion.  It  should  be  noted  that  the  plane  wave,  in  fact,  contains 
all  of  the  functions  from  [22],  p.  767, 

00  00 

(3.25)  expfik  . >(}  = Att  E i'"  Y*^  (o  ' ,<j>' ) . 


m = 0 n = -m 


Here  (0',0')  are  the  spherical  polar  angles  of  jc. 


If  the  source  f or  the  radiated  field  u has  a convergent  expansion 
in  spherical  harmonics,  then  the  convergence  in  (3.22)  is  assured. 
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3.4 


Uniqueness  in  the  Inverse  Source  Problem 


It  has  been  shown  that  the  inverse  source  problem  has 
nonunique  solutions.  The  only  way  to  obtain  a unique  solution  is  to 
provide  additional  information  about  the  source.  Here,  examples  of 
this  type  will  be  presented. 

Example  1:  The  point  source 

Here,  it  is  assumed  a priori  that 

(3.26)  F(x,  t)  = 6(x)  G(t) 

and  that  the  observed  data  is  consistent  with  this  assumption  (i.e.,  an- 
gularly independent).  Although  the  general  theory  will  provide  the 
solution  in  this  case,  it  is  certainly  easier  to  simply  observe  that 

(3.27)  u()^,  w)  = g(w)  exp{iwx/c}/(4TTx) 

and  determine  g uniquely  (hence,  F(x,  t)  uniquely)  from  observations 
in  one  direction. 

A more  interesting  example  is  the  following: 

Example  2:  The  implusive  source 
Here, 

(3.28)  F(x,  t)  = 'S(t)  G(x), 
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The  function  f is  determined  from  observations  by  (3.23),  for  all 
choices  of  the  direction  of  All  values  of  k are  determined  from 
frequency  information.^  Thus  the  function  g(j<)  (and  hence  G {^)) 
is  known.  Because  of  the  a priori  knowledge  (3.28)  about  F(>^,  t), 
this  function  is  now  determined  uniquely. 

This  problem  of  making  the  source  unique  by  additional 
constraint  is  not  well  understood  at  all.  The  authors  believe  that 
a better  understanding  of  the  relation  of  uniqueness  to  side  constraints 
would  have  important  application  to  source  synthesis  and  antenna  design. 


The  problems  of  aperture  limiting  will  not  be  addressed  here. 
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